args <- commandArgs(TRUE)
setwd(args[1])
gene<-read.table("temp_input",sep="\t",head =T)
anno<-read.table("/data1/bsi/BORA_processing/devel/eqtl/parallelize_genotypingsurvival/result_ref_file",sep="\t")
geno<-read.table("/data1/bsi/BORA_processing/devel/eqtl/parallelize_genotypingsurvival/input_nohead.tped",sep=" ")
rownames(geno)<-geno[,2]
rownames(anno)<-as.vector(anno[,1])
gene1<-merge(anno,gene,by="row.names",all.x=F)
i<-1
samplesize<-20
num<-as.integer((length(gene1[1,])-4)/samplesize)
num<-20
predgene<-matrix(nrow=nrow(gene1), ncol=ncol(gene1)-5)
rownames(predgene)<-gene1[,1]
#predgene<-c()
while(i<nrow(gene1)+1)
#while(i<3)
{
        #if(gene1[i,3] == "+"){l1 = gene1[i,4]-30000;l2=gene1[i,4]+20000}else{l1 = gene1[i,4]-20000;l2=gene1[i,4]+30000}
        l1<-gene1[i,4]-20000
	l2<-gene1[i,5]+20000
	l3<-gene1[i,3]
	index<-c()
	index<-which(geno[,4]>=l1 & geno[,4]<l2 & geno[,1] == l3)
	if(length(index) > 0)
	{
        	geno1<-geno[index,]
        	eqtlgene<-as.matrix(gene1[i,6:ncol(gene1)])
        	eqtlgeno<-as.matrix(geno1[5:ncol(geno1)])
		d<-1
		while(d<num+1)
		{			
			s<-d;
			sub<-c()
			main<-c()
			if(d > 1)
			{
				main<-seq(1,d-1,1)
			}
			nu<-1
			while(s<ncol(eqtlgene)+1)
			{
				sub[nu]<-s
				sequ<-c()		
				if(s+1 <= ncol(eqtlgene))
				{	
					sequ<-seq(s+1,s+num-1,1)
					if(s+num-1 > ncol(eqtlgene))
					{
						sequ<-seq(s+1,ncol(eqtlgene),1)
					}
				}
				main<-c(main,sequ)
				nu<-nu+1	
				s<-s+num
			}
			geneeqtl<-t(as.matrix(eqtlgene[,main]))
			geneeqtlsub<-t(as.matrix(eqtlgene[,sub]))
			if(length(index)==1)
			{
                		genomaineqtl<-t(as.matrix(eqtlgeno[,main]))
                		genosubeqtl<-t(as.matrix(eqtlgeno[,sub]))
			}
			else
			{
				genomaineqtl<-eqtlgeno[,main]
                		genosubeqtl<-eqtlgeno[,sub]
			}
			pred1<-c()
			rms<-c()
			u<-1
		
			while(u<nrow(genomaineqtl)+1)
			{
				d1 <- data.frame(geneeqtl[1,],genomaineqtl[u,])
                		colnames(d1)<-c("gene","genotype")
                		glm1 <- glm(gene ~ genotype, data=d1,family=gaussian(link="log"))
                		vec1<-as.data.frame(as.vector(genosubeqtl[u,]))
               	 		colnames(vec1)<-c("genotype")
                		pred <- predict(glm1, data.frame(cbind(genotype=vec1)), type="response")
				pred1<-rbind(pred1,pred)
				rms[u]= 0;
				s<-1
				while(s<length(pred)+1)
				{
					rms[u] = rms[u]+((pred[s]-geneeqtlsub[1,s])*(pred[s]-geneeqtlsub[1,s]));	
					s<-s+1
				}
				u<-u+1	
			}
			nn <-which(rms == min(rms))
			if(length(nn) > 1)
			{
			 	nn<-nn[1]
			}
			t<-1
                	while(t<length(sub)+1)
                	{
                		predgene[i,sub[t]]<-pred1[nn,t]
                		#        #print(pred[i,])
                      		t<-t+1
                	}

			d<-d+1
		}
	}
	print(i)
	i<-i+1
}
colnames(predgene)<-colnames(gene)
write.table(predgene,"predgene_bs_eqtl.txt",quote=FALSE,sep="\t",col.names=TRUE,row.names=TRUE)
